x=randint(3,4,5);dy=rand(3,4);
gamma=rand(1,4);
beta=zeros(1,4);
eps=0.0001;
[xhat,mu,gammax,xmu,ivar,sqrtvar,var,eps,y]...
                =batchnorm_forward(x,gamma,beta,eps);
[dx,dgamma,dbeta]=batchnorm_backward_try...
                (x,y,dy,xhat,mu,gamma,gammax,xmu,ivar,sqrtvar,var,eps);

X=[];dY=[];Y=[];gammaX=[];Xhat=[]
for k=1:4
    X(:,:,k)=x(:,k);
    dY(:,:,k)=dy(:,k);
    gammaX(:,:,k)=gammax(:,k);
    Xhat(:,:,k)=xhat(:,k);
end
G=gamma';
B=beta';
Ybn=vl_nnbnorm(X, G, B);
[dX, dG, dB] = vl_nnbnorm(X, G, B, dY);
[dx(:),dX(:)];
% [y(:),Ybn(:)]


Bn=1;m=Bn;Ki=4;
pre=ivar;
dX1=zeros(size(X));
for ko=1:Ki
    for b=1:Bn
        dy=dY(:,:,ko,b);
        tm1=dY(:,:,ko,:);
        m1=mean(mean(mean(tm1)));
        tm2=dY(:,:,ko,:).*Xhat(:,:,ko,:);
        m2=mean(mean(mean(tm2)));
        M2=m2*gammaX(:,:,ko,b);
        dX1(:,:,ko,b)=pre(ko)*(dy-m1-M2);
    end
end

[dX1(:) dX(:)];
